Mechanism of the Verwey transition in magnetite 
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By combining ab initio results for the electronic structure and phonon spectrum with the group 
theory, we establish the origin of the Verwey transition in FeaCU. Two primary order parameters 
with X$ and As symmetries are identified. They induce the phase transformation from the high- 
temperature cubic to the low-temperature monoclinic structure. The on-site Coulomb interaction 
U between 3d electrons at Fe ions plays a crucial role in this transition — it amplifies the coupling 
of phonons to conduction electrons and thus opens a gap at the Fermi energy. 
Published in Phys. Rev. Lett. 97, 156402 (2006). 

PACS numbers: 71.30.+h, 71.38.-k, 64.70.Kb, 75.50.Gg 



h 

CO 

c3 



i 

C 

O 
o 



> 

m 

m 
o 

o 



o 
o 



X 



The discovery of the remarkable discontinuous drop of 
the electrical conductivity by two orders of magnitude in 
magnetite (Fe^C^) below the Verwey transition (VT) at 
Ty = 122 K [l| triggered intensive studies of its micro- 
scopic origin which have been continued for more than 
six decades. Verwey suggested explaining it by a metal- 
insulator transition due to the reduction of electron mo- 
bility caused by ordering of Fe 3+ and Fe 2+ ions below 
Ty. In spite of great experimental effort, however, the 
existence of ionic-like charge ordering at T < Ty could 
not be confirmed, and the origin of the VT in magnetite 
remains still puzzling 0]. 

Magnetite is a ferrimagnetic spinel with anomalously 
high critical temperature T c ~ 860 K. Hence, it is viewed 
as an ideal candidate for room temperature spintronic 
applications. It crystalizes in the inverse spinel cubic 
structure, FdSm, with two types of Fe sites: the tetra- 
hedral A sites and the octahedral B ones, see Ref. 0. 
Below T c , Fea04 is magnetically ordered with antipar- 
allel moments at A and B sites. Within the commonly 
used ionic picture A sites are occupied by Fe 3+ ions, and 
B sites are occupied randomly by Fe 3+ and Fe 2+ ions at 
room temperature. The Verwey model assumes a purely 
electronic mechanism of the VT, leading below Ty to or- 
dering of Fe 2+ and Fe 3+ ions in B-chains, along [110] and 
[110] directions, respectively. Anderson £| argued that 
the charge ordering in a form of a charge density wave, is 
stabilized by interionic electrostatic energy, where each 
J5-tetrahedron consists of two Fe 2+ and two Fe 3+ ions. 
This picture was studied in the Hubbard-like model with 
interatomic d — d Coulomb interactions [j| , but was not 
confirmed by nuclear magnetic resonance and Mossbauer 
measurements While recent resonant soft X-ray scat- 
tering suggests charge-orbital ordering within the oxygen 
2p states Q and a cooperative Jahn- Teller effect 0, the 
microscopic mechanism of the VT remains controversial. 

There are also other arguments against simple ionic 
mechanism of the VT. Firstly, the replacement of oxygen 
O 16 by O 18 isotope increases Ty by a few degrees in- 
dicating that the transition cannot be purely electronic. 



Secondly, signals in diffuse neutron scattering observed 
above Ty at kA = (0, 0, A), halfway between T and X 
points with kx = (0, 0, 1), and equivalent points of the 
reciprocal lattice [lfj (in units of ^C, where a is the cubic 
lattice constant), change into Bragg peaks below Ty. The 
intensity of critical scattering was succesfully calculated 
on the basis of the transverse acoustic (TA) phonon mode 
with the A5 symmetry. Further neutron studies discov- 
ered diffuse scattering with large intensities close to T and 
X points, and the dominant component of A3 symm etry 
|ll| . Other observations from Raman [J^i X-ray [lflj 
and nuclear inelastic scattering measurements [l^ sug- 
gest that phonons play an essential role in the VT. But 
phonons alone cannot explain the metal-insulator nature 
of the transition either. Actually, if they alone were re- 
sponsible, a phonon soft mode would have been found for 
the cubic phase, while the present ab initio calculations 
do not imply such a soft mode behavior. 

The purpose of this Letter is to resolve the above con- 
troversies and to demonstrate a cooperative scenario of 
the VT. Following the pioneering work of Ihle and Lorenz 
, who considered weak intersite Coulomb interactions 
lfij . we argue that a combination of local Coulomb in- 
teractions between 3d electrons and the electron-phonon 
(EP) coupling is the key feature responsible for the ob- 
served VT. Before presenting the evidence we recall that 
the monoclinic structure at low temperature [J3> space 
group P2/c, with unit cell close to aj\[2 x a/V2 x 2a, 
was observed using the high resolution neutron and syn- 
chrotron X-ray diffraction Q . Having P2/c symmetry of 
the low temperature phase, we get a consistent descrip- 
tion of the VT, with reflections arising below Ty at recip- 
rocal lattice point kA- Reflections at kjf and equivalent 
points are related to the unit cell twice shorter, namely 
a/ V2 x a/v2 x a. Under typical circumstances one would 
select A5 symmetry as a primary order parameter (OP), 
while A3 as a secondary one. The symmetry analysis, 
however, prohibits the A3 symmetry to be a secondary 
OP of A 5 . 

The electronic structures of magnetite in the cubic and 
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monoclinic magnetically ordered phases |3j] have been ob- 
tained using ab initio methods. In the cubic phase, the 
ground state is metallic with the minority spin t 2g Fe(B) 
states at the Fermi energy [l^. The calculations for the 
monoclinic P2/c structure stable below Ty, performed 
using the LDA+U method jyjj, revealed the insulating 
state with a small gap of 0.18 eV, the value being remark- 
ably close to the experimental one 0.15 eV j^- The gap 
opening is accompanied by rather subtle charge-orbital 
ordering, with the charge modulation amplitude in 3d 
Fe(-B) states of order O.le, exactly the same as the one 
derived from diffraction analysis |3j], and below the sen- 
sibility of other experimental techniques . 

In order to combine the observed and computed prop- 
erties of magnetite to a consistent picture of the VT, we 
have performed a symmetry analysis based on the group 
theory. For this study, we used two computer codes: the 
COPL and Isotropy We recall that the high- 
and low-symmetry space groups for cubic and monoclinic 
phases Q are Fd3m and P2/c, respectively. The first ob- 
servation is that there is no single irreducible represen- 
tation (IR) of Fd3m which lowers the space group Fd3m 
directly to P2/c (see Tab. I). This means that the VT has 
to involve at least two OPs. Further symmetry analysis 
shows that two IRs, A 5 and A 3 , acting simulataneously, 
reduce the space group Fd3m to P2/c. Indeed, there are 
two symmetry reduction relationships: Fd3m — > A5 — > 
Pbcm(4) and Fd3m — » A3 — > Pmna(2), where the in- 
crease of the unit cells is indicated in brackets. Common 
symmetry elements of Pbcm(4) and Pmna(2) form the 
space group P2/c(4). (We have verified that the sym- 
metry elements are correctly located and oriented.) As 
a highly nontrivial result one finds that two independent 
IRs, A 5 and A 3 , give two primary OPs of the VT. More- 
over, the Tj = IR could also be involved in the VT. 



TABLE I: List of OPs (IR) from the parent space group Fd3m 
(No=227) and basis (1, 0, 0), (0, 1, 0), (0, 0, 1) to the mono- 
clinic phase P2/c (No=13, unique axis 6, choice 1), basis 
(i,-i,0),(|,|,0), (0,0,2) and origin (±,0, \) relative to the 
original face center cubic lattice. Size is the ratio of primitive 
low-symmetry to high-symmetry unit cell volumes. Specific 
relationship between components of the order parameters are 
not shown. 
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FIG. 1: (color online) Low-frequency phonon dispersion rela- 
tions for cubic phase of Fe3C>4 calculated with (7 = 4 eV and 
J = 0.8 eV. The squares show the results of neutron scattering 
experiment [2*sjl . Two primary OPs are related to A5 and X3 
phonons marked by circles. The high symmetry points from 
left to right are: L = (§, 1, 1), T = (0,0,0), X = (0,0,1), 

k = r = (1,1,1) 



Its symmetry reduction relationship is Fd3m — ► — ► 
Imma(l). Symmetry elements of Imma(l) will not re- 
duce further the symmetry of the space group P2/c(4). 
Thus, is a secondary OP. It represents a shear and 
contributes to the VT, as suggested by the observed soft- 
ening of the C44 elastic constant 23] . 

The group theory tells us the order of coupling terms 
between the OPs. One finds that linear A5 8X3, linear- 
bilinear A5 <8>A3 (&X3 and A5 ® A5 ^1X3 terms are forbid- 
den by symmetry. Therefore, the fourth order coupling 
A5 ® A5 ® A3 ® A3 is the lowest one allowed by symme- 
try. The phase transition with two OPs has to be of the 
first order. The coupling between the secondary OP 
and A3 is described by the linear-bilinear (g> A3 <g> A3 
term. There are other IRs from T, X and A reciprocal 
lattice points, which could become active at the VT as 
secondary OPs. Although they might further diminish 
the ground state energy, they cannot lower the crystal 
symmetry. 

The crystal structure of magnetite was optimized 
within the generalized gradient approximation (GGA) 
with on-site 3d electron interactions described by the 
Coulomb element U and the Hund's exchange J [24j , the 
so-called spin-polarized GGA+J7 approach. The calcu- 
lations were performed with the VASP code on tw° 
supercells: ax ax a and approximately a/V2xa/\/2x 2a, 
for the cubic and monoclinic structures, respectively, 
with 56 atoms each. We included six valence electrons for 
oxygen (2s 2 2p 4 ) and eight for iron (3<i 7 4s 1 ) represented 
by plane waves with energy cut-off 520 eV. The wave 
functions in the core region were obtained by the full- 
potential projector augmented- wave method |26j . The 
summation over the Brillouin zone was performed on the 
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FIG. 2: (color online) Total energies Etot per supercell of 56 
atoms for the structures: cubic Fd3m, cubic with superim- 
posed phonon modes of A5 at Ica (Fd3m+A5), and X3 at 
kx (Fd3m+X3) for increasing mode amplitude Q, and for 
the ground state monoclinic P2/c. Parameters as in Fig. 1. 

4x4x4 and 4x4x2 k-point grids for Fd3m and P2/c 
symmetries, respective ly. For the local interactions we 
choose the parameters 27] : U = 4 eV and J = 0.8 eV. 

Phonons were calculated only for the cubic phase, us- 
ing the direct method implemented in Phonon code |28j . 
The Hellmann-Feynman forces were obtained for six inde- 
pendent displacements: two for each nonequivalent atom 
in positive and negative directions (with the amplitude 
0.02 A). Using the respective force constants, the dy- 
namical matrix was constructed and diagonalized. The 
selected lxlxl supercell provides exact phonon frequen- 
cies at the r and X points. The frequencies away from 
these points carry only negligible errors since the force 
constants decrease more than two orders of magnitude 
within the supercell. 

The phonon dispersion curves, along the high- 
symmetry directions of the reciprocal space, are classi- 
fied according to their IRs (Fig. 1). The longitudinal 
acoustic and TA modes in [001] direction correspond to 
one-dimensional Ai and two-dimensional A5 represen- 
tations, respectively. As already mentioned, all phonon 
frequencies are real, and soft modes are absent in the 
cubic phase. The longitudinal acoustic modes are in a 
very good agreement with the neutron data measured at 
room temperature |29j |. The discrepancies for the trans- 
verse phonons may result from the EP coupling, which 
effectively lowers phonon frequencies at low temperatures 
(as observed by diffuse scattering). 

The OPs A5 and X3 could be (each one separately) a 
linear combination of the electronic and phononic com- 
ponents. In order to investigate the EP coupling for these 
modes, we computed the effect of applied lattice deforma- 
tion, generated according to phonon polarization vectors, 
on the total energy and electron density of states. Among 
all the investigated modes with the A symmetry at kA, 
only the TA (A5) mode shows a significant coupling to 
electrons (presented below). At the X point, we found 
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FIG. 3: (color online) Electron minority density of states for 
the structures of Fig. 2: cubic Fd3m, cubic with superim- 
posed phonon modes of A5 at kA (Fd3m+As) and X3 at kx 
(Fd3m+X3) symmetry with normal coordinate Q = 3, and 
monoclinic P2/c structures. Solid (dashed) lines correspond 
to U = 4 eV and J = 0.8 eV (U = J = 0), respectively. 

a strong EP interaction for the lowest transverse optic 
mode with the A3 symmetry. In contrast, the TA (A4) 
and higher optic phonons couple very weakly to electrons. 

When the cubic crystal is distorted by either X 3 or A 5 
phonons, the total energy E tot of the system decreases, 
indicating that these modes participate in the structural 
transition (Fig. 2). Note that if electrons were not in- 
volved, the ground state energies of the optimized struc- 
tures distorted by phonons would have increased. The 
lowest energy was reached for the P2/c structure, with 
the optimized lattice constants and atomic positions close 
to these of Ref. , what confirms that monoclinic sym- 
metry is stable at T = K. 

The crystal distortion which decreases the energy E to t 
is directly connected with the opening of the gap at Ep 
in the electronic density of states (Fig. 3). For two dis- 
torted structures, and for cubic Fd3m and monoclinic 
P2/c ones, we compare the electronic structures obtained 
for U = and U — 4 eV. In agreement with the previous 
calculations 0, only down spin t 2g Fc(B) states con- 
tribute at Ep , and the cubic phase is metallic indepen- 
dently of U. For the A5 deformation, magnetite remains 
a metal for U — 0, while for U = 4 eV one observes a sig- 
nificant reduction of the spectral weight above Ep, which 
does not yet lead to a gap. This indicates however a con- 
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siderable enhancement of the EP coupling for the increas- 
ing Coulomb interaction U. In the case of the X$ mode 
and for {7 = 4 cV, the EP coupling is even much stronger 
— it opens a gap and triggers the metal-insulator tran- 
sition. In contrast, a metallic phase is found for U = 0. 
Finally, for the monoclinic P2 /c symmetry, the electronic 
state is insulating (metallic) for U = 4 eV (U — 0), in 
agreement with other calculations [lfSj . 

Altogether, the above results reveal a crucial role 
played by electron correlations in the VT. The strong 
Coulomb interaction U reduces the mobility of electrons 
in t 2g states, increasing their tendency towards localiza- 
tion. So modified electronic density responds to lattice 
deformation, leading to the electronic instability at the 
Fermi surface and to the gap opening. We emphasize 
that the present mechanism is novel and does not bene- 
fit from the Peierls-like distortions. Such a cooperative 
mechanism including strong electron correlations and the 
EP coupling was studied before [sO] > but the dimer-bond 
formation seems not to be supported experimentally 

As argued by Wright et al. [3j, the lattice distor- 
tion and charge density following from diffraction analy- 
sis have predominantly character of X modulation, with 
additional A modulation. This is in perfect agreement 
with our result, which shows a strong coupling of the X3 
phonon mode to the electronic density. In this mode, 
the Fe(-B) and O atoms are displaced along the diagonal 
[110] and [110] directions, creating a polar deformation 
of B site octahedra. This lattice distortion couples to 
charge fluctuations on Fe(_B) and O ions, inducing the 
diffuse scattering observed much above Ty. The observa- 
tion that lattice deformation survives locally in the cubic 
phase 0] indicates the existence of the precursor short- 
range order above Ty. It is further supported by photoe- 
mission spectroscopy, consistent with the reduction of the 
single particle gap to ~ 0.1 eV, not closing completely at 
T > Ty [2(J. The long-range order, with larger gap, sets 
in at Ty due to crystal-structural transformation, as was 
clearly demonstrated in recent high-pressure diffraction 
experiment \?>\\ . The observed charge disproportionation 
results from the modified d-p hybridization due to the 
EP coupling. Finally, a gap in the magnon spectrum at 
kA point below Ty 32] indicates a large spin-phonon in- 
teraction and confirms that Ica becomes a point on the 
Brillouin zone surface. 

Summarizing, we have shown that the VT is driven by 
two primary OPs A5 and Y3, which are both character- 
ized by strong linear EP coupling, amplified by Coulomb 
interactions. While the cubic phase would remain metal- 
lic without this coupling, charge fluctuations of two pri- 
mary OPs induced by phonons simultaneously support 
local lattice distortions and open a pseudogap at the 
Fermi energy. As a result, the condensation of both OPs 
leads to the monoclinic distortion below the VT, a gap 
develops and the conductivity is lowered. At the struc- 
tural transition, the charge modulation with a tiny ampli- 



tude is stabilized, being a manifestation of the common 
electronic and lattice origin of this transition. 

Note added. After submitting this paper new resonant 
X-ray scattering data were published [33 which reveal 
fractional charge ordering below the VT. They are fully 
compatible with the present explanation of the VT. 

We acknowledge partial support by the European 
Community under FP6 contract No. NMP4-CT-2003- 
001516 (DYNASYNC). 
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